*******Boot-strap***********




*Set root data directory
local rootdir
cd "`rootdir'"



****Prep Rural-urban****


local list grains veg

local agg regioncode sector
local agg2 regioncode
local baseround=66
local base x


use round`baseround'_edit, clear
drop *fru*


gen lnexp2=log(totexp-xfood)
drop *cloth* *bev* *intox* *light* *sugar* *proc* *meat* *oil* *milk* *food*
drop h*grains_* h*pulses_* h*veg_* 
cap: drop freemeal*
cap: drop z
cap: drop mkt*
drop if totexp==0 | totexp==.


merge 1:1 round sector subround region subsample segment substratum fsu hhno using round`baseround'_hh, keepusing(mpce hhsize religion dmratioa dfratioa dmratios dfratios nic1 htype hgroup mult headed)
drop _merge

merge m:m region round using mult`baseround', keepusing(regioncode)
drop _merge

drop region


drop if mpce>10000 | mpce<10
drop if hhsize>14 | hhsize==. | hhsize==0


merge m:1 round `agg' using pindex_rur_revised
keep if _merge==3
drop _merge




foreach j in grains pulses veg{
gen rx`j'=x`j'/pindex`j'
}

sort `agg'

foreach i in `list'{
levelsof n`i'tot, local(nmax)
egen hhbase`i'=rowmax(x`i'_1-x`i'_`nmax')

gen basep`i'=qpindex`i'

gen cinbase`i'=0
replace cinbase`i'=1 if n`i'>0 & x`i'>0
}
drop qgrains* qpulse* qveg*

summ base* *inbase*, detail




gen lnexp=log(mpce*hhsize)





***household demographic controls
gen lhhsize=log(hhsize)
gen rm=(dmratioa+dmratios)/hhsize
gen rf=(dfratioa+dfratios)/hhsize

***agriculture dummy
gen agg=0
replace agg=1 if nic1=="0"
replace agg=1 if sector==1 & (htype==1 | htype==4)

*hgroup
gen scst=1
replace scst=0 if hgroup==4 | hgroup==9
*religion
gen nonhindu=0
replace nonhindu=1 if religion==1

**pds
cap: gen pds_qrice=0
cap: gen pds_qwheat=0
replace pds_qrice=0 if pds_qrice==.
replace pds_qwheat=0 if pds_qwheat==.

gen pdsqgrains=pds_qrice+pds_qwheat


replace pds_xrice=0 if pds_xrice==.
replace pds_xwheat=0 if pds_xwheat==.
gen pdsxgrains=pds_xrice+pds_xwheat
gen pdssharegrains=pdsxgrains/xgrains

cap: gen pdsanygrains=0
cap: replace pdsanygrains=1 if pdsxgrains>0

gen pdsqveg=0

**home share
foreach i in `list'{
gen hshare`i'=hx`i'/x`i'
}


***village fixed effects...
egen cluster=group(round sector subround regioncode fsu)
bysort cluster: gen count=_N
drop if count<4 | count>10


****For eps,F regression
foreach i in `list'{
gen logx`i'=log(x`i')
gen logn`i'=log(n`i')
}




gen all=1

save temp_ru, replace



***Prep over time***

local list grains veg 
local agg regioncode round
local agg2 regioncode
local baseround=38
local otherround=66
local base x


use round`baseround'_edit, clear
append using round`otherround'_edit


gen lnexp2=log(totexp-xfood)
drop *cloth* *bev* *intox* *light* *sugar* *proc* *meat* *oil* *milk* *food*
drop h*grains_* h*pulses_* h*veg_* h*fru_*
cap: drop freemeal*
drop z mkt*
drop if totexp==0 | totexp==.


merge 1:1 round sector subround region subsample segment substratum fsu hhno using round`baseround'_hh, keepusing(mpce hhsize religion dmratioa dfratioa dmratios dfratios nic1 htype hgroup headed)
drop _merge
merge 1:1 round sector subround region subsample segment substratum fsu hhno using round`otherround'_hh, keepusing(mpce hhsize religion dmratioa dfratioa dmratios dfratios nic1 htype hgroup headed) update
drop _merge

merge m:m region round using mult`baseround', keepusing(regioncode)
drop _merge
merge m:m region using mult`otherround', keepusing(regioncode) update
drop _merge

drop region


drop if mpce>10000 | mpce<10
drop if hhsize>14 | hhsize==. | hhsize==0



merge m:1 `agg' using pindex_r_revised
keep if _merge==3
drop _merge


drop *fru*





sort `agg'

foreach i in `list'{
levelsof n`i'tot, local(nmax)
egen hhbase`i'=rowmax(x`i'_1-x`i'_`nmax')
gen basep`i'=qpindex`i'
gen cinbase`i'=0
replace cinbase`i'=1 if n`i'>0 & x`i'>0
}
drop qgrains* qpulse* qveg* 



gen lnexp=log(mpce*hhsize)





***household demographic controls
gen lhhsize=log(hhsize)
gen rm=(dmratioa+dmratios)/hhsize
gen rf=(dfratioa+dfratios)/hhsize

***agriculture dummy
gen agg=0
replace agg=1 if nic1=="0"
replace agg=1 if sector==1 & (htype==1 | htype==4)

*hgroup
gen scst=1
replace scst=0 if hgroup==4 | hgroup==9
*religion
gen nonhindu=0
replace nonhindu=1 if religion==1

**pds
replace pds_qrice=0 if pds_qrice==.
replace pds_qwheat=0 if pds_qwheat==.

gen pdsqgrains=pds_qrice+pds_qwheat

replace pds_xrice=0 if pds_xrice==.
replace pds_xwheat=0 if pds_xwheat==.
gen pdsxgrains=pds_xrice+pds_xwheat
gen pdssharegrains=pdsxgrains/xgrains

cap: gen pdsanygrains=0
cap: replace pdsanygrains=1 if pdsxgrains>0


gen pdsqveg=0

**home share
foreach i in `list'{
gen hshare`i'=hx`i'/x`i'
}


***village fixed effects...
egen cluster=group(round sector subround regioncode fsu)
bysort cluster: gen count=_N
drop if count<4 | count>10




****For eps,F regression
foreach i in `list'{
gen logx`i'=log(x`i')
gen logn`i'=log(n`i')
}


gen all=1
save temp_time, replace








**********************************
************Bootstrap code*********
**********************************
set seed 250782

clear
set obs 1000
gen b=_n
gen meangrains=2.16
gen sdgrains=0.128
gen sigmagrains=rnormal(meangrains,sdgrains)

gen meanveg=1.988
gen sdveg=0.194
gen sigmaveg=rnormal(meanveg,sdveg)

gen all=1

save sigmas, replace

summ sigma*





***Rural-urban***

clear
save bsresults_ru, replace emptyok


local list grains veg

local agg regioncode sector
local agg2 regioncode
local baseround=66
local base x



forvalues b=1(1)1000{


use temp_ru, clear

gen b=`b'
merge m:1 b using sigmas
keep if _merge==3

egen z=group(`agg')
egen maxz=max(z)
local maxz=maxz[1]



bsample, strata(z)


gen rep=`b'


foreach i in `list'{
cap: drop logpn`i'_hh
gen logpn`i'_hh=log((x`i'/hhbase`i')^(1/(sigma`i'-1)))
}



xtset cluster

local pntype hh
local samp cinbase

foreach i in `list'{
local controls lhhsize rm rf agg scst nonhindu pds*`i' hshare`i'


gen ec`i'_a=.
gen int`i'_a=.
gen psi`i'_a=.


forvalues j=1(1)`maxz'{
**OLS
capture: reg logn`i' logx`i' `controls' if z==`j' & `samp'`i'==1
capture: replace ec`i'_a=_b[logx`i'] if z==`j'
capture: replace int`i'_a=_b[_cons] if z==`j'


capture: reg logpn`i'_`pntype' logn`i'  `controls' if z==`j' & `samp'`i'==1
capture: replace psi`i'_a=_b[logn`i'] if z==`j'
}
}

foreach i in `list'{
replace x`i'=. if x`i'==0
replace n`i'=. if n`i'==0

bysort `agg': egen x`i'_p50=median(x`i')
bysort `agg': egen x`i'_p90=pctile(x`i'), p(90)
bysort `agg': egen x`i'_p10=pctile(x`i'), p(10)

bysort `agg': egen n`i'_p50=median(n`i')
bysort `agg': egen n`i'_mean=mean(n`i')
bysort `agg': egen n`i'_p90=pctile(n`i'), p(90)
bysort `agg': egen n`i'_p10=pctile(n`i'), p(10)
}

duplicates drop `agg', force
*what else to keep? prices and such for decomposition? median expenditures and expenditure percentiles?
keep `agg' ec* int* psi* qpindex* basep* x*_p* n*_p* n*_mean rep

drop if xgrains_p50==0
sort regioncode sector

foreach i in `list'{
gen p`i'=qpindex`i'

foreach j in a{


replace psi`i'_`j'=. if psi`i'_`j'>10

gen eps`i'_`j'=psi`i'_`j'+(1/ec`i'_`j')

*which price to use -- base?

gen logF`i'_`j'=(log(psi`i'_`j'/eps`i'_`j'))-((eps`i'_`j'-psi`i'_`j')*int`i'_`j')-log(p`i')
gen F`i'_`j'=exp(logF`i'_`j')

*for COL later
gen eta`i'_`j'=psi`i'_`j'/eps`i'_`j'
gen const`i'_`j'= ( (eta`i'_`j'^(eta`i'_`j'/(1-eta`i'_`j'))) - (eta`i'_`j'^(1/(1-eta`i'_`j'))) )^(eta`i'_`j'-1)

forvalues k=10(40)90{
gen u0`i'_`k'_`j'=(x`i'_p`k'/(p`i' * (F`i'_`j'^(eta`i'_`j')) * const`i'_`j') )^(1/(1-eta`i'_`j'))
}


*predicted n
gen n`i'_p50_`j'= ((x`i'_p50*psi`i'_`j')/(eps`i'_`j'*F`i'_`j'*p`i'))^(1/(eps`i'_`j'-psi`i'_`j'))

}
}


gen base=1 if sector==1
foreach i in `list'{

by regioncode: egen bp`i'=min(p`i'*base)
by regioncode: egen bx`i'_p50=min(x`i'_p50*base)
by regioncode: egen bx`i'_p10=min(x`i'_p10*base)
by regioncode: egen bx`i'_p90=min(x`i'_p90*base)

by regioncode: egen bn`i'_p50=min(n`i'_p50*base)
by regioncode: egen bn`i'_mean=min(n`i'_mean*base)
by regioncode: egen bn`i'_p10=min(n`i'_p10*base)
by regioncode: egen bn`i'_p90=min(n`i'_p90*base)

gen relx`i'_p50=log(x`i'_p50/bx`i'_p50)-log(p`i'/bp`i')

foreach j in a{

foreach k in psi`i' eps`i' F`i' n`i'_p50 u0`i'_10 u0`i'_50 u0`i'_90{
by regioncode: egen b`k'_`j'=min(`k'_`j'*base)
}

gen beta`i'_`j'=bpsi`i'_`j'/eps`i'_`j'
gen bconst`i'_`j'=( (beta`i'_`j'^(beta`i'_`j'/(1-beta`i'_`j'))) - (beta`i'_`j'^(1/(1-beta`i'_`j'))) )^(beta`i'_`j'-1)


gen reln`i'_p50_`j'=log(n`i'_p50_`j'/bn`i'_p50_`j')
gen relx`i'_p50_`j'=relx`i'_p50*(1/(eps`i'_`j'-psi`i'_`j'))
gen relpsi`i'_p50_`j'=log(psi`i'_`j'/bpsi`i'_`j')
replace relpsi`i'_p50_`j'=relpsi`i'_p50_`j'*(1/(eps`i'_`j'-psi`i'_`j'))
gen logx0`i'_`j'=log((bx`i'_p50*bpsi`i'_`j')/(bF`i'_`j'*bp`i'*beps`i'_`j'))
replace relpsi`i'_p50_`j'=relpsi`i'_p50_`j'+(logx0`i'_`j'  * ( (1/(eps`i'_`j'-psi`i'_`j')) -  (1/(eps`i'_`j'-bpsi`i'_`j')) ) )
gen relvc`i'_p50_`j'=log(beps`i'_`j'/eps`i'_`j')+log(bF`i'_`j'/F`i'_`j')
replace relvc`i'_p50_`j'=relvc`i'_p50_`j'*(1/(eps`i'_`j'-psi`i'_`j'))
replace relvc`i'_p50_`j'=relvc`i'_p50_`j'+( logx0`i'_`j' * ( (1/(eps`i'_`j'-bpsi`i'_`j')) -  (1/(beps`i'_`j'-bpsi`i'_`j')) ) )


foreach k in 10 50 90{
gen newx`i'_`k'_`j'=(bu0`i'_`k'_`j'^(1-eta`i'_`j')) * p`i' * (F`i'_`j'^(eta`i'_`j')) * const`i'_`j'
gen ratio`i'_`k'_`j'=newx`i'_`k'_`j'/bx`i'_p`k'
gen vargain`i'_`k'_`j'=1-(ratio`i'_`k'_`j'*(bp`i'/p`i'))
gen bnewx`i'_`k'_`j'=(bu0`i'_`k'_`j'^(1-beta`i'_`j')) * p`i' * (F`i'_`j'^(beta`i'_`j')) * bconst`i'_`j'
gen bratio`i'_`k'_`j'=bnewx`i'_`k'_`j'/bx`i'_p`k'
gen bvargain`i'_`k'_`j'=1-(bratio`i'_`k'_`j'*(bp`i'/p`i'))
}

gen vargdiff`i'_`j'=vargain`i'_90_`j'-vargain`i'_10_`j'
gen bvargdiff`i'_`j'=bvargain`i'_90_`j'-bvargain`i'_10_`j'

gen ces`i'_`j'=1-((n`i'_p50^(-psi`i'_`j'))/(bn`i'_p50^(-bpsi`i'_`j')))
gen bces`i'_`j'=1-((n`i'_p50^(-bpsi`i'_`j'))/(bn`i'_p50^(-bpsi`i'_`j')))

gen mces`i'_`j'=1-((n`i'_mean^(-psi`i'_`j'))/(bn`i'_mean^(-bpsi`i'_`j')))
gen mbces`i'_`j'=1-((n`i'_mean^(-bpsi`i'_`j'))/(bn`i'_mean^(-bpsi`i'_`j')))

}

}




**variance decomposition
foreach i in `list'{
foreach j in a{


gen vall`i'_`j'=0
foreach k in x psi vc{
bysort sector: egen vrel`k'`i'_p50_`j'=sd(rel`k'`i'_p50_`j')
replace vall`i'_`j'=vall`i'_`j'+(vrel`k'`i'_p50_`j'^2)
}
foreach k in x psi vc{
replace vrel`k'`i'_p50_`j'=(vrel`k'`i'_p50_`j'^2)/vall`i'_`j'
}

}
}

keep if sector==2
keep regioncode *ces* *vargain* *vargdiff* regioncode reln* relpsi* relx* relvc* rep
append using bsresults_ru
save bsresults_ru, replace
}






*****Boot-strap for over time


clear
save bsresults_time, replace emptyok



local list grains veg 
local agg regioncode round
local agg2 regioncode
local baseround=38
local otherround=66
local base x



forvalues b=1(1)1000{


use temp_time, clear

gen b=`b'
merge m:1 b using sigmas
keep if _merge==3

egen z=group(`agg')
egen maxz=max(z)
local maxz=maxz[1]



bsample, strata(z)

gen rep=`b'



foreach i in `list'{
cap: drop logpn`i'_hh
gen logpn`i'_hh=log((x`i'/hhbase`i')^(1/(sigma`i'-1)))
}



xtset cluster

local pntype hh
local samp cinbase

foreach i in `list'{
local controls lhhsize rm rf agg scst nonhindu pds*`i' hshare`i'

*a - ols 

gen ec`i'_a=.

gen int`i'_a=.


gen psi`i'_a=.
.

forvalues j=1(1)`maxz'{
**OLS
capture: reg logn`i' logx`i' `controls' if z==`j' & `samp'`i'==1
capture: replace ec`i'_a=_b[logx`i'] if z==`j'
capture: replace int`i'_a=_b[_cons] if z==`j'


capture: reg logpn`i'_`pntype' logn`i' `controls' if z==`j' & `samp'`i'==1
capture: replace psi`i'_a=_b[logn`i'] if z==`j'
}

}


foreach i in `list'{
replace x`i'=. if x`i'==0
replace n`i'=. if n`i'==0
bysort `agg': egen x`i'_p50=median(x`i')
bysort `agg': egen x`i'_p90=pctile(x`i'), p(90)
bysort `agg': egen x`i'_p10=pctile(x`i'), p(10)

bysort `agg': egen n`i'_mean=mean(n`i')
bysort `agg': egen n`i'_p50=median(n`i')
bysort `agg': egen n`i'_p90=pctile(n`i'), p(90)
bysort `agg': egen n`i'_p10=pctile(n`i'), p(10)
}

duplicates drop `agg', force
*what else to keep? prices and such for decomposition? median expenditures and expenditure percentiles?
keep `agg' ec* int* psi* qpindex* basep* x*_p* n*_p* n*_mean round rep


drop if xgrains_p50==0

sort regioncode round

foreach i in `list'{

gen p`i'=basep`i'

foreach j in a{


replace psi`i'_`j'=. if psi`i'_`j'>10

gen eps`i'_`j'=psi`i'_`j'+(1/ec`i'_`j')


gen logF`i'_`j'=(log(psi`i'_`j'/eps`i'_`j'))-((eps`i'_`j'-psi`i'_`j')*int`i'_`j')-log(p`i')
gen F`i'_`j'=exp(logF`i'_`j')

*for COL later
gen eta`i'_`j'=psi`i'_`j'/eps`i'_`j'
gen const`i'_`j'= ( (eta`i'_`j'^(eta`i'_`j'/(1-eta`i'_`j'))) - (eta`i'_`j'^(1/(1-eta`i'_`j'))) )^(eta`i'_`j'-1)

forvalues k=10(40)90{
gen u0`i'_`k'_`j'=(x`i'_p`k'/(p`i' * (F`i'_`j'^(eta`i'_`j')) * const`i'_`j') )^(1/(1-eta`i'_`j'))
}

*predicted n
gen n`i'_p50_`j'= ((x`i'_p50*psi`i'_`j')/(eps`i'_`j'*F`i'_`j'*p`i'))^(1/(eps`i'_`j'-psi`i'_`j'))

}
}




gen base=1 if round==38
foreach i in `list'{

by regioncode: egen bp`i'=min(p`i'*base)
by regioncode: egen bx`i'_p50=min(x`i'_p50*base)
by regioncode: egen bx`i'_p10=min(x`i'_p10*base)
by regioncode: egen bx`i'_p90=min(x`i'_p90*base)

by regioncode: egen bn`i'_p50=min(n`i'_p50*base)
by regioncode: egen bn`i'_p10=min(n`i'_p10*base)
by regioncode: egen bn`i'_p90=min(n`i'_p90*base)

by regioncode: egen bn`i'_mean=min(n`i'_mean)

gen relx`i'_p50=log(x`i'_p50/bx`i'_p50)-log(p`i'/bp`i')

foreach j in a{

foreach k in psi`i' eps`i' F`i' n`i'_p50 u0`i'_10 u0`i'_50 u0`i'_90{
by regioncode: egen b`k'_`j'=min(`k'_`j'*base)
}

gen beta`i'_`j'=bpsi`i'_`j'/eps`i'_`j'
gen bconst`i'_`j'=( (beta`i'_`j'^(beta`i'_`j'/(1-beta`i'_`j'))) - (beta`i'_`j'^(1/(1-beta`i'_`j'))) )^(beta`i'_`j'-1)

gen reln`i'_p50_`j'=log(n`i'_p50_`j'/bn`i'_p50_`j')
gen relx`i'_p50_`j'=relx`i'_p50*(1/(eps`i'_`j'-psi`i'_`j'))
gen relpsi`i'_p50_`j'=log(psi`i'_`j'/bpsi`i'_`j')
replace relpsi`i'_p50_`j'=relpsi`i'_p50_`j'*(1/(eps`i'_`j'-psi`i'_`j'))
gen logx0`i'_`j'=log((bx`i'_p50*bpsi`i'_`j')/(bF`i'_`j'*bp`i'*beps`i'_`j'))
replace relpsi`i'_p50_`j'=relpsi`i'_p50_`j'+(logx0`i'_`j'  * ( (1/(eps`i'_`j'-psi`i'_`j')) -  (1/(eps`i'_`j'-bpsi`i'_`j')) ) )
gen relvc`i'_p50_`j'=log(beps`i'_`j'/eps`i'_`j')+log(bF`i'_`j'/F`i'_`j')
replace relvc`i'_p50_`j'=relvc`i'_p50_`j'*(1/(eps`i'_`j'-psi`i'_`j'))

replace relvc`i'_p50_`j'=relvc`i'_p50_`j'+( logx0`i'_`j' * ( (1/(eps`i'_`j'-bpsi`i'_`j')) -  (1/(beps`i'_`j'-bpsi`i'_`j')) ) )


foreach k in 10 50 90{
gen newx`i'_`k'_`j'=(bu0`i'_`k'_`j'^(1-eta`i'_`j')) * p`i' * (F`i'_`j'^(eta`i'_`j')) * const`i'_`j'
gen ratio`i'_`k'_`j'=newx`i'_`k'_`j'/bx`i'_p`k'
gen vargain`i'_`k'_`j'=1-(ratio`i'_`k'_`j'*(bp`i'/p`i'))

gen bnewx`i'_`k'_`j'=(bu0`i'_`k'_`j'^(1-beta`i'_`j')) * p`i' * (F`i'_`j'^(beta`i'_`j')) * bconst`i'_`j'
gen bratio`i'_`k'_`j'=bnewx`i'_`k'_`j'/bx`i'_p`k'
gen bvargain`i'_`k'_`j'=1-(bratio`i'_`k'_`j'*(bp`i'/p`i'))

}

gen vargdiff`i'_`j'=vargain`i'_90_`j'-vargain`i'_10_`j'
gen bvargdiff`i'_`j'=bvargain`i'_90_`j'-bvargain`i'_10_`j'


gen ces`i'_`j'=1- ((n`i'_p50^(-psi`i'_`j'))/(bn`i'_p50^(-bpsi`i'_`j')))
gen bces`i'_`j'=1- ((n`i'_p50^(-bpsi`i'_`j'))/(bn`i'_p50^(-bpsi`i'_`j')))

gen mces`i'_`j'=1- ((n`i'_mean^(-psi`i'_`j'))/(bn`i'_mean^(-bpsi`i'_`j')))
gen mbces`i'_`j'=1- ((n`i'_mean^(-bpsi`i'_`j'))/(bn`i'_mean^(-bpsi`i'_`j')))


}

}



keep if round==66
keep regioncode *ces* *vargain* *vargdiff* regioncode reln* relpsi* relx* relvc* rep
append using bsresults_time
save bsresults_time, replace
}












****Analysis
use bsresults_ru, clear





foreach k in grains veg{
foreach j in reln`k'_p50_a relx`k'_p50_a relpsi`k'_p50_a relvc`k'_p50_a bvargain`k'_50_a vargain`k'_50_a bvargdiff`k'_a vargdiff`k'_a{
bysort regioncode: egen mean_`j'=mean(`j')
bysort regioncode: egen sd_`j'=sd(`j')
}
}

keep mean_* sd_* regioncode
duplicates drop regioncode, force


set scheme s2mono 
local k grains
gen tstat`k'=mean_bvargain`k'_50_a/sd_bvargain`k'_50_a
gen atstat`k'=mean_vargain`k'_50_a/sd_vargain`k'_50_a
twoway (kdensity tstat`k', xline(1.96) xline(-1.96) xlabel(-1.96 "-1.96" 1.96 "1.96") title("Grains: t-statistics by region") xtitle(t-statistic) ytitle(Density) legend(order(1 "Welfare gain (variety cost only)" 2 "Welfare gain (all)"))) (kdensity atstat`k', xline(1.96) xline(-1.96))
graph save `k'_tstat, replace
twoway (kdensity mean_bvargain`k'_50_a, title("Grains: welfare point estimates by region") xline(0) xtitle("Welfare gain (compensating variation, % of category expenditure)") ytitle(Density) legend(order(1 "Welfare gain (variety cost only)" 2 "Welfare gain (all)"))) (kdensity mean_vargain`k'_50_a)
graph save `k'_estimate, replace
graph combine `k'_estimate.gph `k'_tstat.gph, rows(2)
graph export `k'_regions_ru.png, replace

replace mean_vargainveg_50_a=. if mean_vargainveg_50_a<-2
replace mean_bvargainveg_50_a=. if mean_bvargainveg_50_a<-2

local k veg
gen tstat`k'=mean_bvargain`k'_50_a/sd_bvargain`k'_50_a
gen atstat`k'=mean_vargain`k'_50_a/sd_vargain`k'_50_a
twoway (kdensity tstat`k', xline(1.96) xline(-1.96) xlabel(-1.96 "-1.96" 1.96 "1.96") title("Veg: t-statistics by region") xtitle(t-statistic) ytitle(Density) legend(order(1 "Welfare gain (variety cost only)" 2 "Welfare gain (all)"))) (kdensity atstat`k', xline(1.96) xline(-1.96))
graph save `k'_tstat, replace
twoway (kdensity mean_bvargain`k'_50_a, title("Veg: welfare point estimates by region") xline(0) xtitle("Welfare gain (compensating variation, % of category expenditure)") ytitle(Density) legend(order(1 "Welfare gain (variety cost only)" 2 "Welfare gain  (all)"))) (kdensity mean_vargain`k'_50_a)
graph save `k'_estimate, replace
graph combine `k'_estimate.gph `k'_tstat.gph, rows(2)
graph export `k'_regions_ru.png, replace


*****Figure A.13****
graph combine grains_estimate.gph grains_tstat.gph veg_estimate.gph veg_tstat.gph, altshrink
graph export regions_ru.pdf, replace









use bsresults_time, clear

foreach k in grains veg{
foreach j in reln`k'_p50_a relx`k'_p50_a relpsi`k'_p50_a relvc`k'_p50_a bvargain`k'_50_a vargain`k'_50_a bvargdiff`k'_a vargdiff`k'_a{
bysort regioncode: egen mean_`j'=mean(`j')
bysort regioncode: egen sd_`j'=sd(`j')
}
}

keep mean_* sd_* regioncode
duplicates drop regioncode, force





set scheme s2mono 

replace mean_vargaingrains_50_a=. if mean_vargaingrains_50_a<-2
replace mean_bvargaingrains_50_a=. if mean_bvargaingrains_50_a<-2
local k grains
gen tstat`k'=mean_bvargain`k'_50_a/sd_bvargain`k'_50_a
gen atstat`k'=mean_vargain`k'_50_a/sd_vargain`k'_50_a
summ mean_bvargain`k'_50_a mean_vargain`k'_50_a, detail



twoway (kdensity tstat`k', xline(1.96) xline(-1.96) xlabel(-1.96 "-1.96" 1.96 "1.96") title("Grains: t-statistics by region") xtitle(t-statistic) ytitle(Density) legend(order(1 "Welfare gain (variety cost only)" 2 "Welfare gain (all)"))) (kdensity atstat`k', xline(1.96) xline(-1.96))
graph save `k'_tstat, replace
twoway (kdensity mean_bvargain`k'_50_a, title("Grains: welfare point estimates by region") xline(0) xtitle("Welfare gain (compensating variation, % of category expenditure)") ytitle(Density) legend(order(1 "Welfare gain (variety cost only)" 2 "Welfare gain (all)"))) (kdensity mean_vargain`k'_50_a)
graph save `k'_estimate, replace

local k grains
graph combine `k'_estimate.gph `k'_tstat.gph, rows(2)
graph export `k'_regions_time.png, replace



replace mean_vargainveg_50_a=. if mean_vargainveg_50_a<-2
replace mean_bvargainveg_50_a=. if mean_bvargainveg_50_a<-2

local k veg
gen tstat`k'=mean_bvargain`k'_50_a/sd_bvargain`k'_50_a
gen atstat`k'=mean_vargain`k'_50_a/sd_vargain`k'_50_a
summ mean_bvargain`k'_50_a mean_vargain`k'_50_a, detail

twoway (kdensity tstat`k', xline(1.96) xline(-1.96) xlabel(-1.96 "-1.96" 1.96 "1.96") title("Veg: t-statistics by region") xtitle(t-statistic) ytitle(Density) legend(order(1 "Welfare gain (variety cost only)" 2 "Welfare gain (all)"))) (kdensity atstat`k', xline(1.96) xline(-1.96))
graph save `k'_tstat, replace
twoway (kdensity mean_bvargain`k'_50_a, title("Veg: welfare point estimates by region") xline(0) xtitle("Welfare gain (compensating variation, % of category expenditure)") ytitle(Density) legend(order(1 "Welfare gain (variety cost only)" 2 "Welfare gain (all)"))) (kdensity mean_vargain`k'_50_a)
graph save `k'_estimate, replace
graph combine `k'_estimate.gph `k'_tstat.gph, rows(2)
graph export `k'_regions_time.png, replace

****Figure A.12*****
graph combine grains_estimate.gph grains_tstat.gph veg_estimate.gph veg_tstat.gph, altshrink
graph export regions_time.pdf, replace 






*****Revised code for generating table with full region by region boot-strapped means and standard errors
**9 columns, with s.e. underneath


**Over time
use bsresults_time, clear


foreach k in grains veg{
foreach j in reln`k'_p50_a relx`k'_p50_a relpsi`k'_p50_a relvc`k'_p50_a bvargain`k'_50_a vargain`k'_50_a bvargdiff`k'_a{
bysort regioncode: egen sd_`j'=sd(`j')
}
}
keep sd_* regioncode
duplicates drop regioncode, force

*merge in point estimates
merge 1:1 regioncode using ecv_results3866, keepusing(reln*_p50_a relx*_p50_a relpsi*_p50_a relvc*_p50_a bvargain*_50_a vargain*_50_a bvargdiff*_a)
drop _merge

*merge in ces results
merge 1:1 regioncode using cesreg_3866, keepusing(lambdagrains lambdaveg)
replace lambdagrains=-lambdagrains
replace lambdaveg=-lambdaveg
drop _merge
merge 1:1 regioncode using cesvill_3866, keepusing(medianlambdagrains medianlambdaveg)

replace medianlambdagrains=-medianlambdagrains
replace medianlambdaveg=-medianlambdaveg
drop _merge

*merge in region names
merge 1:1 regioncode using region_data, keepusing(stateregionname)
drop _merge
drop if regioncode==.
replace stateregionname="Gujarat Plains Southern" if stateregionname==""

*order and sort to create s.e.

foreach k in grains veg{
preserve
*local k grains

keep stateregionname regioncode *`k'*

gsort - reln`k'_p50_a
gen rank=_n

expand 2
bysort regioncode: gen counter=_n



foreach j in reln`k' relx`k' relpsi`k' relvc`k'{
gen `j'=.
replace `j'=`j'_p50_a if counter==1
replace `j'=sd_`j'_p50_a if counter==2
drop `j'_p50_a sd_`j'_p50_a
}

foreach j in bvargain`k' vargain`k'{
gen `j'=.
replace `j'=`j'_50_a if counter==1
replace `j'=sd_`j'_50_a if counter==2
drop `j'_50_a sd_`j'_50_a
}

foreach j in bvargdiff`k'{
gen `j'=.
replace `j'=`j'_a if counter==1
replace `j'=sd_`j'_a if counter==2
drop `j'_a sd_`j'_a
}

foreach j in lambda`k' medianlambda`k'{
replace `j'=. if counter==2
}

drop if reln`k'==.
bysort regioncode: gen obs=_N
drop if obs==1

sort rank counter

replace stateregionname="" if counter==2
order stateregionname reln`k' relx`k' relpsi`k' relvc`k' bvargain`k' vargain`k' bvargdiff`k' lambda`k' medianlambda`k'
browse stateregionname reln`k' relx`k' relpsi`k' relvc`k' bvargain`k' vargain`k' bvargdiff`k' lambda`k' medianlambda`k'
keep  stateregionname reln`k' relx`k' relpsi`k' relvc`k' bvargain`k' vargain`k' bvargdiff`k' lambda`k' medianlambda`k'
***Appendix: Table A.10 and A.11
save regions_`k'_time, replace
restore
}



**Rural urban
use bsresults_ru, clear

foreach k in grains veg{
foreach j in reln`k'_p50_a relx`k'_p50_a relpsi`k'_p50_a relvc`k'_p50_a bvargain`k'_50_a vargain`k'_50_a bvargdiff`k'_a{
bysort regioncode: egen sd_`j'=sd(`j')
}
}
keep sd_* regioncode
duplicates drop regioncode, force

*merge in point estimates
merge 1:1 regioncode using ecv_results66ru, keepusing(reln*_p50_a relx*_p50_a relpsi*_p50_a relvc*_p50_a bvargain*_50_a vargain*_50_a bvargdiff*_a)
drop _merge

*merge in ces results
merge 1:1 regioncode using cesreg_66ru, keepusing(lambdagrains lambdaveg)
replace lambdagrains=-lambdagrains
replace lambdaveg=-lambdaveg
drop _merge
merge 1:1 regioncode using cesvill_66ru, keepusing(medianlambdagrains medianlambdaveg)

replace medianlambdagrains=-medianlambdagrains
replace medianlambdaveg=-medianlambdaveg
drop _merge

*merge in region names
merge 1:1 regioncode using region_data, keepusing(stateregionname)
drop _merge
drop if regioncode==.
replace stateregionname="Gujarat Plains Southern" if stateregionname==""

*order and sort to create s.e.

foreach k in grains veg{
preserve
*local k grains

keep stateregionname regioncode *`k'*

gsort - reln`k'_p50_a
gen rank=_n

expand 2
bysort regioncode: gen counter=_n



foreach j in reln`k' relx`k' relpsi`k' relvc`k'{
gen `j'=.
replace `j'=`j'_p50_a if counter==1
replace `j'=sd_`j'_p50_a if counter==2
drop `j'_p50_a sd_`j'_p50_a
}

foreach j in bvargain`k' vargain`k'{
gen `j'=.
replace `j'=`j'_50_a if counter==1
replace `j'=sd_`j'_50_a if counter==2
drop `j'_50_a sd_`j'_50_a
}

foreach j in bvargdiff`k'{
gen `j'=.
replace `j'=`j'_a if counter==1
replace `j'=sd_`j'_a if counter==2
drop `j'_a sd_`j'_a
}

foreach j in lambda`k' medianlambda`k'{
replace `j'=. if counter==2
}

drop if reln`k'==.
bysort regioncode: gen obs=_N
drop if obs==1
replace stateregionname="" if counter==2

sort rank counter

***Appendix: Table A.12 and A.13
order stateregionname reln`k' relx`k' relpsi`k' relvc`k' bvargain`k' vargain`k' bvargdiff`k' lambda`k' medianlambda`k'
browse stateregionname reln`k' relx`k' relpsi`k' relvc`k' bvargain`k' vargain`k' bvargdiff`k' lambda`k' medianlambda`k'
keep  stateregionname reln`k' relx`k' relpsi`k' relvc`k' bvargain`k' vargain`k' bvargdiff`k' lambda`k' medianlambda`k'
save regions_`k'_ru, replace
restore
}
